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Abstract 


A finite-difference, three dimensional incompressible Navier-Stokes formulation 
to calculate the flow through turbopump components is utilized. The solution method 
is based on the pseudocompressibility approach and uses an implicit-upwind differenc- 
ing scheme together with the Gauss-Seidel line relaxation method. Both steady and 
unsteady flow calculations can be perfomed using the current algorithm. In this work, 
the equations are solved in steadily rotating reference frames by using the steady-state 
formulation in order to simulate the flow through a turbopump inducer. Eddy viscosity 
is computed by using an algebraic mixing-length turbulence model. Numerical results 
are compared with experimental measurements and a good agreement is found between 
the two. Time-accurate calculations, such as impeller and diffusor interaction, will be 
reported in future work. 


Introduction 


With the advent of supercomputer hardware as well as fast numerical meth- 
ods, computational fluid dynamics (CFD) has become an essential part of aerospace 
research and design. Numerical studies in incompressible flows show good progress 
in parallel with computational studies in compressible flows. For example, the incom- 
pressible flow solver developed by Kwak et al [1] was extensively used for simulating the 
flow through space shuttle main engine power head components. The redesign of the 
space shuttle main engine hot gas manifold, guided by the computations of Chang et 
al. [2], illustrates the usefulness of CFD in aerospace research. Since the incompressible 
Navier-Stokes formulation does not yield the pressure field explicitly from the equation 
of state or through the continuity equation, numerical solution of the equations requires 
special attention in order to satisfy the divergence-free constraint on the velocity field. 
The most widely used methods which use primitive variables are fractional-step and 
pseudocompressibility techniques. In the fractional-step method, the auxiliary veloc- 
ity field is solved by using the momentum equations. Then, a Poisson equation for 
pressure is formed by taking the divergence of the momentum equations and by using 
a divergence-free velocity field constraint. Solving the Poisson equation for pressure 
efficiently in three-dimensional curvilinear coordinates is the most important feature of 
the fractional step method. 3 One way to avoid the numerical difficulty originated by 
the elliptic nature of the problem is to use a pseudocompressibility method. With the 
pseudocompressibility method, the elliptic-parabolic type equations are transformed 
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into hyperbolic-parabolic type equations. Well established solution algorithms devel- 
oped for compressible flows can be utilized to solve the resulting equations. 

Steger and Kutler 4 employed an alternating direction implicit scheme into 
Chorin’s 5 pseudocompressibility method. This formulation was extended to three- 
dimensional generalized coordinates by Kwak. 1 Recently, a three-dimensional incom- 
pressible Navier-Stokes solver (INS3D-LU-SGS) using a lower-upper symmetric-Gauss- 
Seidel algorithm was developed by Yoon and Kwak. 6 This algorithm is used to calculate 
the inducer flow of the Space Shuttle Main Engine turbopump in order to demonstrate 
the performance of the numerical method. 7 Another effort is performed in Ref. 8 by 
using upwind differencing and Gauss- Seidel line relaxation scheme in order to have a 
robust and fast converging scheme (INS3D-UP). A time accurate formulation of this 
algorithm is implemented for incompressible flows through artificial heart devices with 
moving boundaries. 8 ’ 10 In the present study, the steady-state formulation is used in 
steadily rotating reference frames in order to develop a CFD procedure for simulating 
the flow through turbopump components of a liquid rocket engine. 


Computed Results 


The flowfield through a turbopump inducer is solved as a benchmark problem in 
order to validate the CFD procedure for turbomachinery applications. In this section 
results obtained for the Rocketdyne inducer shown in Fig. 1 are presented. The inducer 
geometry was developed and experimentally studied by the Rocketdyne Division of 
Rockwell International. The design flow is 2236 GPM with a design speed of 3600 
RPM. In the computational study, tip-leakage effects are included with a tip clearance 
of 0.008 inches. The problem was nondimensionalized with a tip diameter of 6.0 inches 
and the average inflow velocity of 28.3 ft/sec. The Reynolds number for this calculation 
was 191,800. The upstream section of the inducer was taken as a two tip-diameter- 
long straight channel, as shown in Fig. 1. The bull-nose of the inducer was treated 
as a rotating wall and the cavity section was neglected. However, this region can be 
included by using an additional zone. An H-H grid topology with dimensions of 187 
x 27 x 35 was used. A partial view of the surface grid is shown in Fig. 2. An H-type 
surface grid was generated for each surface using an elliptic grid generator. The interior 
region of the three-dimensional grid was filled using an algebraic solver coupled with 
an elliptic smoother. In the straight channel, the grid was generated for one-sixth of 
the cross-section of the tube. This grid was extended to the outflow section of the 
inducer between the blades. Periodic boundary conditions were used at the end points 
in the rotational direction. At inflow and outflow boundaries characteristic boundary 
conditions were employed. At the inflow, v and w velocity components were specified as 
zero and the total pressure was specified as constant. Axial velocity and static pressure 
were calculated from the characteristic relation and the total pressure relation. At 
the outflow, static pressure was specified and the velocity components were computed 
from the characteristics propagating from the interior region. The flow was taken at 
rest initially and the inducer was fully rotated impulsively. The solution was considered 
■ converged when the maximum residual dropped at least four orders of magnitude. This 
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was obtained in less than 500 iterations. Computer time required per grid point per 
iteration was about 1.4 x 10 -4 sec. 

Figure 3 illustrates the planes where the experimental measurements were taken 
by Rocketdyne. Axial and tangential velocity components and the flow angle were 
measured in planes A,B,C and D at various circular arcs from the hub to the tip region. 
At each plane, the comparison between experimental measurements and numerical 
results along three of the circular arcs is presented in this paper. A total velocity and 
a flow angle are compared against experimental data. The total velocity has only a 
tangential and an axial velocity components. The radial velocity component was not 
measured in the experiment. 

Figures 4 through 7 show relative total velocities and relative flow angles as 
a function of circumferential angle in degrees in planes A, B, C, and D, respectively. 
The circumferential angle increases from the suction side to the pressure side. The 
dashed lines in these figures represent the experimental data and solid lines represent 
the numerical results. The comparison of computations and experiment is generally 
good all the way from the hub to tip region. The difference between experimental and 
numerical data is about 5-8 % in velocity. In all planes, the hub and tip regions indicate 
the biggest discrepancy. This may be a result of the relatively coarse grid used for the 
boundary layer. In the computational study, the Baldwin-Lomax algebraic turbulence 
model is used to determine the eddy viscosity. The comparison shows that the solution 
algorithm does a god job with an algebraic turbulence model. The implementation of 
the one equation model 11 of Baldwin and Barth is currently underway for the present 
algorithm. The motivation for higher-order turbulence modeling is due to the compar- 
isons obtained in Plane D, in which the wake region is not predicted accurately (Fig. 
7). Another advantage of the one-equation model is that there is no need to define 
a length-scale explicitly. Near the tip clearance region, the difference between experi- 
mental measurements and numerical results is noticably larger than the error in other 
regions. This is due to lack of grid resolution in the tip clearance region. In the grid 
refinement study, the number of grid points in the tip clearance region was increased 
from 4 points to 9 points. In the coarse grid computation, there is one overlapped 
grid point in the rotational direction to ensure periodic boundary conditions. In the 
fine grid, additional 3 zones were added in radial direction. The results with the one 
equation model and the results from the grid refinement study will be published in 
future. 

Figure 8 shows the surface of the inducer colored by nondimensionalized pres- 
sure. The pressure gradient across the blades due to the action of centrifugal force and 
the pressure rise from inflow to outflow are illustrated. This pressure rise along the 
inducer can also be seen in Fig. 9. Velocity vectors are plotted in the meridional plane 
and the vectors are colored by the static pressure. The existing solution procedure can 
be applied to the same configuration under off-design conditions. The massive sepa- 
ration which may block the fuel supply can be detected in the numerical study. This 
is the future research area of the present study which can be used as a pre-design and 
post-design engineering tool in challenging turbomachinery applications. 
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Summary 


An efficient and robust solution procedure is implemented and validated for 
three-dimensional turbopump applications. Numerical simulations of the flow through 
the Rocketdyne inducer have been successfully carried out by using CFD techniques for 
solving viscous incompressible Navier-Stokes equations with the source terms in steadily 
rotating reference frames. The method of artificial compressibility with a higher-order 
accurate upwind differencing and the Gauss-Seidel line relaxation scheme provide fast 
convergence and robustness. Results in the form of relative total velocity and relative 
flow angle in four planes are presented. Numerical results compare fairly well with 
experimental data. 
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Figure 1. Rocket dyne turbopump inducer configuration. 



Figure 2: Surface grid for Rocketdyne turbopump inducer. 
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Figure 4: Comparison of relative total velocity and relative flow angle in Plane A 
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Figure 5: Comparison of relative total velocity and relative flow angle in Plane B. 
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Figure 6: Comparison of relative total velocity and relative flow angle in Plane C. 
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Figure 8: Surface pressure for Rocketdyne inducer. 
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Figure 9: Velocity vectors colored by pressure on the meridional plane of the inducer. 
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Abstract. The extension of computational fluid dynamics techniques to 
artificial heart flow simulation is illustrated. Unsteady Navier-Stokes equa- 
tions written in three-dimensional generalized curvilinear coordinates are 
solved iteratively at each physical time-step until the incompressibility con- 
dition is satisfied. The solotion method is based upon the pseudocompress- 
ibility approach and uses an implicit upwind differencing scheme together 
with the Gauss-Seidel line-relaxation method. The efficiency and robust- 
ness of the time-accurate formulation of the numerical algorithm are tested 
by computing the flow through model geometries. A channel flow with 
a moving indentation is computed and validated with experimental mea- 
surements and other numerical solutions. In order to handle the geomet- 
ric complexity and the moving boundary problems, a zonal method and 
an overlapped grid- embedding scheme are employed, respectively. Steady- 
state solutions for the flow tliroufh a tilting-disk heart valve are compared 
with experimental measurements; good agreement is obtained. The flow 
computation during opening and closing of the valve is carried out to il- 
lustrate the moving-boundary capability. Aided by experimental evidence, 
the flow through an entire Penn State artificial heart model is computed. 


1. Introduction 

With the advent of supercomputer hardware, as well as fast numerical meth- 
ods, researchers in the field of computational fluid dynamics (CFD) tackle more 
complicated problems than ever before. With these new capabilities, CFD has 
become an essential part of aerospace research and design. For example, the in- 
compressible flow solver developed by Kwak et al. (Ref. 1) was extensively used 
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State electric artificial heart, can be well defined (Ref. 9). Therefore, the present 
study is focused on the fluid problem with prescribed body motion. 

Underwood and Mueller obtained the flow characteristics for the Kay-Shiley 
disk-type valve by using the stream function-vorticity formulation (Ref. 10). 
Their results showed agreement with experimental data up to a Reynolds number 
of 600. Idelsohn et al. modeled the flow through the Kay-Shiley caged disk, 
Starr-Edwards caged ball, and Bjork-Shiley tilting-disk valves and compared 
their performance (Ref. 11). Turbulent flow through trileaflet aortic heart valves 
was simulated by Steven sen et al. (Ref. 12). Most numerical studies assumed 
that the flow through the heart valve was two-dimensional. Additionally, the 
valve opening and closing motions were neglected; only the flow through a fixed 
valve position was studied. In reality, the geometry is three-dimensional, and 
the flow through heart valves involves moving boundaries. 

In our present study we propose the development of a computational pro- 
cedure simulating steady and unsteady three-dimensional flows through artifi- 
cial hearts and heart valves with moving boundaries. In the next sections, the 
method of solution is summarized followed by a demonstration of how this CFD 
procedure can be used for probing the flow through artificial heart' devices. 

2. Method of solution 

The flow through artificial heart devices is unsteady, viscous, and incompress- 
ible. In the present study, the non-Newtonian nature of the blood is neglected, 
and the flow is described by the three-dimensional incompressible Navier-Stokes 
equations. Numerical simulation of such flows is a very challenging problem in 
computational fluid dynamics. In addition to the geometric complexity, obtain- 
ing time-dependent solutions of the incompressible Navier-Stokes equations with 
moving boundaries poses many difficult numerical problems. 

Because the pressure and velocity fields are not directly coupled owing to 
the lack of a pressure term in the continuity equation, numerical solution of 
the incompressible Navier-Stokes equations requires special attention in order 
to satisfy the divergence-free constraint on the velocity field. The most widely 
used methods that use primitive variables are fractional-step and pseudocom- 
pressibility techniques. In the fractional-step method, the auxiliary velocity field 
is solved by using the momentum equations. Then, a Poisson equation for pres- 
sure is formed by taking the divergence of the momentum equations and by using 
a divergence-free velocity field constraint. The efficient solution of the Poisson 
equation for pressure in three-dimensional curvilinear coordinates is the most 
important feature of the fractional-step method (Ref. 13). One way to avoid 
the numerical difficulty originated by the elliptic nature of the problem is to use 
a pseudocompressibility method. With the pseudocompressibility method, the 
elliptic-parabolic-type equations are transformed into hyperbolic-parabolic-type 
equations. Well-established solution algorithms developed for compressible flows 
can then be utilized to solve the resulting equations. 
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In the present study, the pseudo compressibility approach is used and the time- 
accuracy is attained by iterating in pseudo-time until the divergence of velocity 
is driven toward zero to within a specified tolerance. Here, the time-derivatives 
in the momentum equations are differenced by using a second-order, three-point, 
backward-difference formula. The numerical method uses a second-order central- 
difference for viscous terms and a higher-order flux- difference splitting for the 
convective terms. The resulting matrix equation is solved iteratively by using a 
nonfactored line-relaxation scheme, which maintains stability and allows a large 
pseudo-time-step to be taken. At each sweep direction, a tridiagonal matrix is 
formed, and off-line terms of the matrix equation are moved to the right-hand 
side of the equation. Details of the governing equations and numerical method 
are given in References 14 and 15. 

One of the biggest difficulties in the simulation of flows in complicated three- 
dimensional configurations is the discretization of the physical domain. The 
problem becomes more severe if one body in the domain of interest moves rela- 
tive to another one, as is seen in the tilting disk valve and the Penn State artificial 
heart geometries. The use of a zonal approach is a practical solution if the grids 
are stationary. For more general applications, a chimera grid-embedding tech- 
nique provides a greater flexibility for the grid motion (Ref. 16). This technique 
is employed for flow computations through a tilting disk valve and the Penn 
State artificial heart model. 


3. Computed results 

One of the goals of this study is to simulate the flow through a realistic model 
of an artificial heart. Since the geometry and the flow physics are complicated, 
the computational procedure is validated by solving several simpler problems 
which characterize the flow in various parts of an artificial heart. As a first 
step, an idealized two-dimensional pump model was chosen to demonstrate the 
capability of the time-accurate formulation under a moving grid-condition. The 
geometry of this model and the computed results are presented in Reference 14. 
Channel flow with an identical wall, the flow through a tilting-disk heart valve, 
and the flow through the Penn State artificial heart model are included in this 
section. 

3.1 Channel flow with a moving indentation. 

Channel flow with an asymmetric oscillating indentation was experimentally 
studied by Pedley and Stephanoff and was numerically simulated by Ralph and 
Pedley (Ref. 18). In the experiment, the channel was was rigid everywhere except 
for the indentation, which is made of a thick rubber membrane. The experiment 
shows that flow is two-dimensional near the midplane, and so the computation 
is done in two dimensions. 

Figure 1 illustrates the instantaneous streamlines plotted at several nondimen- 
sional times for Reynolds number Re = 600 and for Strouhal number St = 0.057. 
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The Reynolds number is based on the channel height a, and on the average ve- 
locity at the entrance of the channel U. The Strouhal number is defined as 

St = af/Uy where / is the oscillation frequency. At the beginning of the cycle, 
the flow downstream of the indentation is parallel to the channel walls. A single 
eddy is formed at the sloping wall of the indentation during the first half of the 
cycle. The streamlines at the core flow are lifted slightly upward as shown in 
Figure la. This is the beginning of the wavy flow patterns of the core flow. 
The formation of a second separated eddy on the opposite wall can be seen in 
Figure lb. In later stages, the double row of eddies along the lower and upper 
wall of the channel is observed. At the end of the cycle, the vortices are swept 
downstream (Fig. lg), and the residual vortices are not strong enough to affect 
the next cycle. In fact, the flow pattern of the first and second cycles are quite 
similar. Consequently, the flow is assumed to be periodic in time, even at the 
first cycle. Another interesting phenomenon observed in the present study, as 
well as in the experimental and other numerical studies, is the eddy-doubling 
which can be seen in Figures lc through le. It occurred in the second, third, 
and fourth vortices from the indentation. In eddy-doubling, a single eddy splits 
into two rotating eddies. 

The time-evolution of the centers of the vortices is compared with experimen- 
tal and other numerical findings. The first four vortices are labeled vortices A, 
B, C, and D, and are shown in Figure 2b. The distance between the indenta- 
tion wall and the center of these vortices is measured from the instantaneous 
streamlines and plotted versus time in Figure 2a. The dashed lines represent 
the present computations. The solid lines denote computational results from 
the fractional-step approach implemented by Rosenfeld (Ref. 19). Dotted lines 
are numerical results of stream-function vorticity formulations from Ralph and 
Redley (Ref. 18). Experimental measurements by Pedley and Stephanoff (Ref. 
17) are represented by square symbols. The agreement between the numerical 
and experimental measurements is fairly good. There is a discrepancy between 
the numerical results and experimental measurements about the location of the 
vortex A. The present results and Rosenfeld’s computations predict the location 
of vortex A to be 0.4 units closer to the indentation than the experimental find- 
ings. However, the locations of the remaining vortices are correctly predicted, if 
the distance. is measured from the center of the vortex A. This underprediction 
of the separation length of vortex A is thought to be caused by an inaccurate 
description of the indentation wall shape. In the present study and in Rosen- 
feld’s study the same grid and the same wall shape were used. Even though 
the solution algorithms of the two computational studies are completely differ- 
ent, the agreement between the numerical results is good. For these reasons, 
only the locations of vortices B, C, and D are compared with the experimental 
measurements. 

3.2 Flow through tilting disk valve. 

The problem of flow through the tilting-disk valve was chosen to develop and 
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validate a procedure that will be used for the valve region of the Penn State 
artificial heart. In the Bjork-Shiley tilting-disk heart valve, the tilting disk is 
placed in front of the sinus region of the human aorta. The aortic root has three 
sinuses about 120° apart. The tilting-disk valve model used in this computation 
is simplified by assuming that the sinus region of the aorta has a circular cross 
section. The cage and struts which hold the free-floating disk inside of the sewing 
ring are not included in the geometry. It is also assumed that the walls have no 
elastic deformation. The channel length is taken to be five aorta diameters long. 
The computational geometry used in these unsteady flow computations is given 
in Figure 3. The disk motion is illustrated by showing three different positions 
of the disk, at angles of 75°, 50°, and 30°, as measured from the centerline of the 
aorta. The tilting disk is allowed to rotate about the horizontal axis, which is 
1/6 of a disk diameter below the center of the disk. Because of this asymmetric 
disk orientation, the flow is three dimensional. 

The chimera grid-embedding technique, which has been successfully used for 
external flow problems, has been employed by using two overlapped grids, as 
shown in Figure 4. Grid 1 occupies the whole region in the aorta from entrance 
to exit, and remains stationary. Grid 2 wraps around the tilting disk, and moves 
with the disk. In the chimera grid-embedding technique, grid points that lie 
within the disk geometry and outside the channel grid are excluded from the 
solution process. These excluded points are called hole points, and the immediate 
neighbors of the hole points are called fringe points. The information is passed 
from one grid to another via fringe and grid boundary points by interpolating the 
dependent variables. Tri-linear interpolation is used in the present computations. 
In order to distinguish the hole and fringe points from regular computational 
points, an IBLANK array is used in the flow solver. For hole, grid boundary, 
and fringe points, IBLANK is set to zero; otherwise it is set to one. In order 
to exclude the hole and grid boundary points from the solution procedure, the 
coefficients of the system of algebraic equations and the right-hand-side terms are 
multiplied by the IBLANK value. If the grid point is a hole, an outer boundary, 
or a fringe point, the value of (1 - IBLANK) is added to the main diagonal of 
the matrix equation. 

Presented here are the results of steady flow with a fixed disk angle and 
unsteady flow with the disk motion in the configuration described above. The 
problems are nondimensionalized by using the entrance diameter as a unit length, 
and the average inflow velocity as the unit velocity. In order to reduce the 
computational effort and memory size, the inflow and outflow boundaries are 
placed a short distance from the region of interest relative to the boundaries used 
in the experimental studies. In addition, the exact shape of the sinus region of 
the aorta used in the experiments is not known. These discrepancies could lead 
to slight differences between the present computations and experimental results. 

Steady-state calculations for the 30° disk orientation have been carried out for 
Reynolds numbers in the range of 2,000 to 6,000, for which experimental data 
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are available. The Reynolds number is based on the diameter and the mean 
velocity at the entrance of the channel. A mixing-length algebraic turbulence 
model, which is used for incompressible flow through the space shuttle main- 
engine turnaround duct (Ref. 2), is utilized. Figure 5 shows the pressure drop 
across the Bjork-Shiley tilting-disk valve at different flow rates of physiological 
interest. The computed and measured axial velocity profiles at 42 mm down- 
stream from the disk are shown in Figure 6. Axial velocity profiles are plotted 
in the horizontal plane through the center of the channel. The numerical results 
are shown with dots, and the experimental results are shown with triangles. The 
numerical results compare favorably with the experimental measurements (Ref. 
3). The largest discrepancy is seen near the walls, where the boundary layer 
is overestimated by the calculation. Figure 7 shows the velocity vectors at five 
longitudinal stations. The flow, which is directed to the upper part of the aorta, 
generates vortices in the sinus region of the aorta and a large separated region 
along the lower wall of the aorta. Since separated-flow and low-flow regions have 
the potential to form thrombi, clotting may occur on the upper sinus region and 
the lower wall of the aorta. Figliola and Mueller also present mean velocity pro- 
files, which show similar flow characteristics to those indicated by computational 
study, at several locations (Ref. 5). They computed the shear stress from the 
measured velocity field and observed that the maximum shear occurs at the top 
wall downstream of the sinus region of the aorta. This is in agreement with the 
velocity plot shown in Figure 7, in which there are large velocity components 
just off the wall in that location. Particle traces in Figure 8 indicate that the 
flow does not separate adjacent to the tilting disk. The tilting disk separates the 
flow into a major flow region, which is along the upper wall of the tube, and a 
minor flow region, which is along the lower wall of the tube. Separation, reverse 
flow, and swirling motion mostly occur in the minor flow region. 

Figure 9 shows vorticity magnitude contours in the surface of the tube, outflow 
surface, and inflow surface of the disk, respectively. It is assumed that maximum 
vorticity magnitudes indicate the regions of high shear. The sewing ring surface 
and the edges of the disk are the regions where the vorticity magnitude is at a 
maximum. The upper wall of the channel also has high vorticity magnitudes. 

Unsteady-flow calculations have been carried out in order to demonstrate and 
analyze the flow during opening and closing of the disks. For the present com- 
putations, one cycle of valve opening and closing requires 70 physical time-steps. 
During each time-step, subiterations are carried out until both the maximum di- 
vergence of velocity and the maximum residual drop below 10" 3 . The computing 
time required for one cycle of the valve opening and closing is approximately 5 
Cray-YMP hours. During the valve opening, inflow velocity is imposed at the 
entrance of the channel The inflow velocity is chosen as a sine function in time. 
The disk rotation is specified as a linear function in time. Since the forces acting 
on the disk are known from the numerical solution, the disk rotation angle can be 
determined. However, the disk rotation angle is limited. For large disk-rotation 
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angle, some information may be lost between the grids when the grid-embedding 
technique is used. In order to prevent the information loss, the maximum allowed 
disk-rotation angle at each physical time-step is taken to be less than 3°. 

Figures 10a through 10c illustrate the velocity vectors on the lateral symmetry 
plane at t/T = 0.128, 0.257, and 0.385, respectively. The velocities are very high 
in the region between the disk and the channel wall, as shown in Figure 10a. 
During the disk opening, two vortices are formed at the upper and lower edges of 
the disk. The flow starts to separate behind the disk and reattaches to the wall as 
shown in Figure 10b. The stagnation region behind the disk moves downstream 
as the disk rotates. Highly skewed velocity profiles are seen downstream from 
the disk, as illustrated in Figure 10c. The growth of the vortices has also been 
observed in the sinus region of the aorta while the flow opens the valve. Along 
the lower wall a separation region is formed. 

3.3 Artificial heart flow. 

The Penn State artificial heart model is composed of a cylindrical chamber 
with two tube extensions (Fig. 11). The inflow (mitral) and outflow (aortic) 
tubes contain concave tilting disks, which open and close to act as valves. In the 
computational model, tilting-disk mitral valve orientation in time was obtained 
from the experimental data provided by Pennsylvania State University. The 
aortic valve orientation in time was approximated to mitral valve orientation 
with a phase difference. The pumping action is provided by a pusher plate 
whose velocity is sinusoidal in time. The diameter of the pusher plate is 7.26 
cm, and it has a stroke length of 2.28 cm. The problem is nondimensionalized 
with the inflow tube diameter, which is 2.54 cm, and a unit velocity of 20 cm/sec. 
In the computational study, the Reynolds number based on the unit length and 
velocity is 900. Initially, the flow was started at rest, and four cycles of the 
pumping action were completed using a Cray-YMP computer at NASA-Ames 
Research Center. One cycle of the pusher plate’s motion required 240 physical 
time-steps. At each time-step, the equations were iterated until the maximum 
divergence of velocity was reduced below 10” 2 . During most of the cycle, 10-20 
subiterations were required (for more detail, see Refs. 14, 15). 

In order to handle the geometric complexity and the moving boundary-pro- 
blems, a zonal method and an overlapped grid-embedding scheme (Ref. 16) are 
employed, respectively. In the zonal method, a complex computational domain is 
divided into several simple subdomains. The overlapped grid-embedding scheme 
allows subdomains to move relative to each other, and provides great flexibility 
when the boundary movement creates large displacements. The computational 
grid for the heart model is shown in Figure 11. Grid 1 is generated for the pusher 
plate and moves with it. Grid 2 occupies the chamber and remains stationary. 
Grids 3 and 5 are for the inflow and outflow tube extensions, respectively. Grid 
points for the tubes and grid points for the chamber are overlapped on three 
common planes. In other words, the grid points for the tubes start three stencils 
inside of the chamber outer boundaries. Zonal boundary conditions are used in 
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the interface boundaries. Grids 4 and 6 wrap around tilting disks, and move with 
the disks. An overlapped grid-embedding scheme is employed between moving 
grids and stationary grids. 

The purpose of the computed results presented next is mainly to demonstrate 
how CFD can be used to understand the flow in the artificial heart, and the com- 
parison with experiment is qualitative. Unsteady particle traces are illustrated 
in Figure 12. The particles are released near the inflow valve at the beginning of 
the fourth cycle. The figure is plotted at nondimensional time t/T = 0.45 into 
the period at which time the pusher plate is close to its lowest position, where 
T denotes the period for the pusher plate's motion. Figure 12 shows that the 
flow creates a strong vortex in the center region of the chamber. The particles 
have a swirling motion against the back wall against the mitral valve opening. 
The flow also separates at the connection region of the chamber and the inflow 
tube. Figure 13 shows the computed velocity vectors on the horizontal midplane 
at nondimensional time 0.375. At that time, the pusher plate is moving down, 
and the mitral valve is open. Figure 13 also indicates the presence of strong 
circulation in the chamber. However, the three-dimensional structure of the flow 
cannot be seen clearly, because the vectors are plotted on a two-dimensional 
plane. 

The strong vortex in the center of the chamber is actually created where the 
chamber and the inflow tube are connected. The vortex moves to the core of 
the chamber in time. Experimental measurements by Baldwin and Tarbel (Ref. 
19) are illustrated in Figure 14. Since the computational study does not include 
the blood sac inside the chamber, the comparison of experimental and compu- 
tational results is qualitative. In addition, the Reynolds number in the experi- 
mental study is 1.7 times larger than the Reynolds number in the computational 
study, because the flow is assumed to be laminar in the present computations. 
The biggest discrepancy between experimental measurements and computational 
results is the location of the vortex core in the chamber. In Figure 13, the vortex 
is off-center in the chamber. In Figure 14, the vortex is located almost in ''the 
center of the chamber. Another difference can be seen in the wake region of the 
mitral valve. Since the Reynolds number in the computational study is lower, 
the wake is not as strong as the wake in Figure 14. The Reynolds number in 
the future computational study will be increased by including the turbulence 
modeling in order to have a quantitative comparison between experimental and 
computational results. 

During the second half of the cycle time, the pusher plate moves upward, and 
the outflow valve is opened. A top view of the computed velocity vectors on 
the horizontal plane at t/T = 0.625 is plotted in Figure 15. Since the inflow 
valve is closed, residual eddies are quite large near the disk. However, they 
are quickly weakened as the pressure builds up inside the chamber. Measured 
velocity vectors at t/T — 0.625 are shown in Figure 16. 

Figure 17 shows additional vortex structures in the vertical plane of the cham- 
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ber at nondimensional time 0.25. That vortex structure causes the swirling mo- 
tion of the particles previously shown in Figure 12. The swirling flow in the 
outflow tube is shown in Figure 18. The velocity vectors in the cross-sectional 
plane of the outflow tube at t/T = 0.75 are plotted. The cross-sectional plane 
is one nondimensional unit downstream of the tilting-disk valve, and the normal 
vector of that cross section is in the positive ^-direction shown in Figure 13. 

4. Summary 

An efficient and robust solution procedure is developed and validated for nu- 
merical simulations for internal flows through artificial heart devices. The so- 
lution procedure for unsteady, incompressible, viscous flow computations has 
been extended with the incorporation of the grid-embedding approach. This has 
been used to simulate the flow through a tilting-disk heart valve and the flow 
through the Penn State artificial heart model. Separated and secondary flow 
regions have been pointed out in the tilting-disk heart valve and artificial heart 
flow simulations. The vortex created in the central portion of the Penn State 
artificial heart provides good wall washing over the entire chamber. The present 
capability of simulating complicated internal flow problems with moving bound- 
aries is demonstrated. The procedure developed in this study is quite general 
and is applicable to various types of artificial heart and valve geometries. It is 
hoped that artificial heart researchers and designers may further benefit from 
the computational ability obtained in the current work. 
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Figure Captions 

Figure 1. Instantaneous streamlines: Re = 600, St — 0.057. (a) t = 0.45. (b) 
t = 0.55. (c) t = 0.65. (d) t = 0.75. (e) t = 0.80. (f) t = 0.90. (g) t = 1.5. 

Figure 2. Comparison of numerical results with experimental measurements 
and other numerical solutions. 

(a) Time-evolution of vortices center. 

(b) Vortices A, B, C, and D. 

Figure 3. Tilting-disk geometry showing valve opening. 

Figure 4. Overlapped grids. 

Figure 5. Pressure drop across tilting-disk valve versus steady-state flow rate. 

Figure 6. Comparison of numerical result and experimental measurements. 

Figure 7. Velocity vectors at several longitudinal stations. 

Figure 8. Particle traces showing minor and major flow regions. 

Figure 9. Vorticity magnitude contours. 

(a) Channel surface. 

(b) Outflow surface of disk. 

(c) Inflow surface of disk. 

Figure 10. Side view of velocity vectors on the lateral symmetry plane showing 
the valve. 

(a) t/T- 0.128. 

(b) t/T= 0.257. 

(c) i/T = 0.385. 

Figure 11. Computational grid for Penn State artificial heart model showing 
zonal and overlapped grid regions. 

Figure 12. Unsteady particle traces at t/T = 0.45 as pusher plate nears 
bottom position. 

Figure 13. Velocity vectors of incoming fluid at t/T = 0.375. Top view in 
horizontal planes through the center of tubes and in horizontal plane 3 mm below 
top surface of chamber. 

Figure 14. Experimental results at t/T = = 0.375 (Ref. 19). 

Figure 15. Velocity vectors of incoming fluid at t/T = 0.625; top view in 
horizontal planes through center of tubes and in horizontal plane 3 mm below 
top surface of the chamber. 

Figure 16. Experimental results at t/T = 0.625 (Ref. 19). 

Figure 17. Velocity vectors of incoming fluid at t/T = 0.25; side view in 
vertical plane through center of inflow tube. 

Figure 18. Velocity vectors in the cross-sectional plane of outflow tube showing 
swirling motion of the fluid: t/T = 0.75. 
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(a) Time evolution of the vortices center. 



(b) Showing the vortides A.B,C, and D. 


Figure 2: Comparison of the numerical results with the experimental measurements and 
other numerical solutions. 
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Figure 3: Tilting disk geometry showing valve opening. 



Figure 4: Perspective view of two overlapped grids. 
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Figure 5: Pressure drop across the tilting disk valve versus steady-state flow rate. 
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Figure 6: Comparison between the numerical result and the experimental measurements. 




Figure 7: Velocity vectors at several longitudinal stations. 



Figure 8: Particle traces showing the minor and the major flow regions 
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Figure 9: Vorticity magnitude contours on the (a) channel surface; (b) outflow surface of 
the disk; (c) inflow surface of the disk. 
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Figure 10: Side view of velocity vectors on the lateral symmetry plane showing the valve 
openning. 



Figure 11: Computational grid for Penn State artificial heart model showing zonal and 
overlapped grid regions. 



Figure 12: Unsteady particle traces at t/T = 0.45 as the pusher plate nears the bottom 
position. 
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Figure 13: Velocity vectors of incoming fluid at t/T = 0.375. Top view in horizontal planes 

through the center of tubes and in horizontal plane 3 mm. below the top surface 
of the chamber. 



Figure 14: Experimental results at t/T = 0 37=; Kv R l j • 
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